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ABSTRACT 

We model a one-dimensional shock-tube using smoothed particle hydrodynamics and 
investigate the consequences of having finite shock-width in numerical simulations 
caused by finite resolution of the codes. We investigate the cooling of gas during 
passage through the shock for three different cooling regimes. 

For a theoretical shock temperature of 10 5 K, the maximum temperature of the 
gas is much reduced. When the ration of the cooling time to shock-crossing time was 
8, we found a reduction of 25 percent in the maximum temperature reached by the 
gas. When the ratio was reduced to 1.2 maximum temperature reached dropped to 
50 percent of the theoretical value. In both cases the cooling time was reduced by a 
factor of 2. 

At lower temperatures, we are especially interested in the production of molecular 
Hydrogen and so we follow the ionization level and H2 abundance across the shock. 
This regime is particularly relevent to simulations of primordial galaxy formation for 
halos in which the virial temperature of the galaxy is sufficiently high to partially re- 
ionize the gas. The effect of in-shock cooling is substantial: the maximum temperature 
the gas reaches compared to the theoretical temperature was found to vary between 
0.15 and 0.81 for the simulations performed, depending upon the strength of the 
shock and the mass resolution. The downstream ionization level is reduced from the 
theoretical level by a factor of between 2.4 and 12.5, and the resulting H2 abundance 
was found to be reduced to a fraction of 0.45 to 0.74 of its theoretical value. 

At temperatures above 10 5 K, radiative shocks are unstable and will oscillate. 
We reproduce these oscillations and find good agreement with the previous work of 
Chevalier and Imamura (1982), and Imamura, Wolff and Durisen (1984). The effect 
of in-shock cooling in such shocks is difficult to quantify, but is undoubtedly present. 

We conclude that extreme caution must be exercised when interpretting the results 
of simulations of galaxy formation. 
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1 INTRODUCTION 

In order for hydrodynamics codes to be able to simulate 
shocks it is necessary that the width of the shock be spread 
over a few mesh cells, or represent a few inter-particle spac- 
ings, depending on whether a grid or particle based code is 
being used. In reality, however, the true width of the shock 
is only a few mean free paths, giving a shock thickness Ax, 
of 



where n is the number density and a is the collisional cross- 
section. A few simple calculations show that the true shock 
thickness is orders of magnitude smaller than the shock 
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thickness obtained in simulations. This allows the possibil- 
ity that gas in simulations can cool artificially while passing 
through a shock, and that temperature-sensitive quantities, 
such as ionization levels and molecular abundances may also 
vary. Any such effects will be purely numerical because of 
the finite shock crossing time in simulations. In the present 
paper will consider a number of possibilities. Firstly, if gas 
is able to cool during a shock the maximum temperature 
the gas will reach after passing through the shock will be 
reduced from that of an adiabatic jump. Where the cool- 
ing function rate decreases with increasing temperature, this 
will then reduce the cooling time of the post-shock gas. Sec- 
ondly, we shall consider the evolution of the ionization level 
through a shock, where the cooling rate is proportional to 
the number density of free electrons. This has important 
consequences for simulations of primordial structure forma- 
tion as the ionization level determines the rate of molecular 
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hydrogen (H2) and Hydrogen-Deuterium (HD)production. 
These two molecules then determine the cooling time of the 
gas down to temperatures less than 100K. Finally, we shall 
also address the question of instability of radiative shocks. 
Imamura et al. (1984) and Chevalier and Imamura (1982) 
found that radiative shocks with a cooling rate per unit vol- 
ume oc p 2 T a (where p is the mass density and T the tem- 
perature) are unstable against perturbations when 01 < 0.4. 
This instability results in periodic oscillations of the position 
of the shock front and of the jump temperature at the shock 
front. We address the nature of these oscillations in the con- 
text of galaxy formation using a cooling function composed 
of the sum of a number of terms, each with a different power 
law. 



2 ANALYTIC PROFILES 

Shocks can be generally classified into three types, depend- 
ing upon the amount of cooling present. In adiabatic shocks 
the shocked gas does not cool significantly during the pe- 
riod of time of interest. Isothermal shocks have the same 
jump conditions at the shock front as an adiabatic shock, 
but the cooling is very rapid so that the downstream gas is 
at the same temperature as the upstream gas. Any cooling 
tail is thus very narrow. In between these two extremes is 
the case where the shocked gas has a cooling time of order 
the time-scale which is of interest in a given situation. 

For an adiabatic shock the conservative form of the fluid 
equations in the absence of sources which change momen- 
tum, mass or energy can be written as 
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where the variables u, x, p, e and P represent velocity, po- 
sition, density, specific energy and pressure respectively, all 
evaluated in the rest frame of the shock. 

If we apply equations ||-^ to gas far upstream and write 
its properties in terms of gas far downstream, we arrive at 
the Rankine-Hugoniot jump conditions: 
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where T is the gas temperature, k b the Boltzmann constant 
p the mean molecular weight and ?tih the mass of a hydro- 
gen atom. The suffixes 1 and 2 refer to the upstream and 
downstream gas respectively and the equation of state of the 
gas is given by the equations: 

P = k b T 
p pniH 

(9) 
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Taking 7 = 5/3 we can write the pressure of the shocked gas 
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For the case of an isothermal shock where the temper- 
ature of the gas far downstream is equal to its temperature 
up stream, the conservation of energy equation is no longer 
valid and should be replaced with the condition Ti = T2. 
We must stress however that the gas will still undergo an 
adiabatic jump at the shock front before cooling back down 
to its upstream temperature. 

For the case of a radiative shock we wish to determine 
the path between the state of the gas at the shock front and 
its state downstream. In order to do this we must add a 
cooling term to the energy flux equation: 
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The cooling term A serves to determine the path between 
the down stream and adiabatic jump conditions, but does 
not change them in any way. Equations |, I and can be 
numerically integrated for any cooling function resulting in 
an analytic profile from the specified down-stream condi- 
tions to the shock front. The shock front is reached when 
the adiabatic jump condition of equation [l(] is satisfied. 



3 SHOCK-TUBE 

3.1 Description of the code 

The shock tube simulations were performed using a 
smoothed particle hydrodynamics code. The code was that 
described by Couchman, Thomas & Pearce (1995) but with 
modified geometry to make the simulation volume cuboidal 
with dimensions of 6 x 6 x 100. The shock front was ini- 
tially set to be at z coordinate 50 with gas flowing in the z 
direction. On what will be the upstream side of this point, 
particles are placed one per unit volume and given unit mass. 
On the downstream side, particles were placed four per unit 
volume. Their mass however is allowed to vary in order for 
the gas to have any density profile that is initially desired. 
The density profile, along with the temperature and veloc- 
ity of the gas, is read in upon the start of a simulation. For 
this geometry the simulation thus starts with 9000 particles. 
The particles on each side of the shock front were allowed 
to relax independently to a glass-like initial state. 

The sides of the simulation box were periodic, with the 
exception of each end, where gas was fed into the box at the 
appropriate rate upstream and removed from the simula- 
tion a sufficient distance downstream so as not to effect any 
cooling tail. A padding region of two mesh cells at each end 
was necessary and gas within this region had its position up- 
dated on each time-step, but did not have its energy, velocity 
or density updated. (The code searches for 32 neighbours, 
thus needing a search radius of approximately 2 units when 
particles are located one per cell. It is this distance that de- 
termines the size of the padding region as any particle closer 
than this distance to the end of the simulation volume would 
experience a non-symmetric force.) All simulations were per- 
formed in the rest frame of the shock and the units the code 
used were dimensionless until a cooling function was chosen. 

The cooling function, and ionization/recombination 
rates were integrated using RK4 from Press et al. (1992), 



© 0000 RAS, MNRAS 000, 000-000 



In-Shock Cooling in Numerical Simulations 3 



this routine being called at the end of each time-step. All 
the rates used were tabulated on the first timestep. 



4 STABLE IN-SHOCK COOLING 
4.1 Theory 

The cooling function for objects with high virial tempera- 
tures (above 10 4 K) collapsing after the Universe has been 
enriched with metals can be written as the sum of three 
terms: a bremsstrahlung term, Ai, a metal line-cooling term, 
A2 and a H-He term, A3. The fit for each term is split into 
two temperature regimes. For T > 10 K we have, where 



A = A/erg cm 3 s 1 : 

Ai = 5.2 x 10~ 28 (T/K)' (12) 

A 2 = 1.7 x lO" 18 (T/K)" a8 Z/Z (13) 

A 3 = 1.4 x 1(T 18 (T/Ky 1 (14) 
For 10 5 K > T > 10 4 K we have: 

Ai=0 (15) 

A 2 = 1.7 x 10" 27 (T/K)Z/Z Q (16) 

A 3 = 1.4 x 1CT 28 (T/K) (17) 



where Z/Z© is the metallicity in terms of solar. 

The above cooling function is a good approximation to 
within a factor of 2 to Raymond & Smith (1977). The na- 
ture of this cooling function is such that the cooling rate 
rises abruptly as one moves from a temperature T < 10 4 K 
to T > 10 4 K. The cooling rate peaks at 10 5 K and then 
decreases, reaching a minimum at 2 x 10 7 K. At this point 
bremsstrahlung cooling starts to become effective and the 
cooling rate begins to slowly rise as the temperature in- 
creases further. For gas being shocked from 10 4 K to a tem- 
perature of order 2 x 10 7 K the cooling function is of a suit- 
able nature to produce a significant reduction in the cooling 
time of shocked gas, if in-shock cooling is present. The cool- 
ing rate being most rapid at lower temperatures and then 
decreasing as the temperature rises above 10 5 K means that 
a modest reduction in the maximum temperature reached 
could make a significant difference to the cooling time as 
it is the hottest gas that takes the longest to cool. Unfortu- 
nately a declining cooling function makes the shock unstable 
and in this section we restrict ourselves to studying shocks 
which do not reach the declining part of the cooling function, 
and which thus settle down into a steady state. 

The nature of such shocks can be specified by the ra- 
tio of the cooling time to the shock crossing time. If the 
cooling time is very long, and the ratio is thus large, the ef- 
fect of in-shock will be negligible. Conversely, if the cooling 
time is very rapid, in-shock cooling will tend to drive the 
shock towards an isothermal transition. The cooling time it- 
self will be determined by the physical conditions that are 
being simulated and will be determined by the density and 
temperature of the gas. The shock crossing time will depend 
upon the mass resolution of the code and will vary as the 
cube root of the particle mass. 
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Figure 1. Steady state temperature profile demonstrating in- 
shock cooling. The solid line shows the analytic profile that would 
be expected in the absence of in-shock cooling. The ratio of the 
cooling-time to shock crossing time is of order 8. 



4.2 Methodology 

By specifying the downstream conditions of the shock, equa- 
tions ^, |^ and [ll] were numerically integrated for A of section 
4.1, and a profile of the cooling tail was obtained. The in- 
tegration was stopped at the shock front, when equation |ici| 
was satisfied. This profile was then read directly into the 
shock-tube to create the initial conditions, with the shock 
front being placed at the boundary between the low and 
high mass particles. Particles within the region of the cool- 
ing tail had their temperature, velocity and mass adjusted 
to recreate the analytic shock profile. The simulations were 
then allowed to run until a steady state-state was reached. 

The upstream gas was started with a temperature 
of 5.95 x 10 4 K and the analytic jump temperature was 
1.84 x 10 5 K. The density of gas far downstream had a den- 
sity twelve times that of the upstream gas. Two simulations 
were ran, one with a ratio of the shock crossing time to the 
cooling time of 8, and the other with a ratio of 1.2. 

4.3 Results 

Figures ^ and ^ show the temperature profiles as a function 
of position for two shocks with ratios of 8 and 1.2 respec- 
tively. Despite the upper part of the theoretical shock being 
in the unstable regime the simulations were found to settle 
down into a steady state. This was aided by in-shock cooling 
which reduces the maximum temperature reached. The solid 
lines in show the analytic temperature profiles. We can see 
that the maximum temperature obtained is approximately 
25 percent less than the analytic jump temperature in Fig- 
ure |l] , and 50 percent less for the simulation of Figure g 
with a particle mass 512 times greater. 

The effect that these reductions in jump temperature 
have on the cooling tails can be seen to be somewhat er- 
ratic. In both cases the cooling time of the gas decreases 
by a factor of approximately 2, but the spatial extent of the 
cooling tail depends upon the resolution. The reason for this 
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Figure 2. Steady state temperature profile demonstrating in- 
shock cooling. The solid line shows the analytic profile that would 
be expected in the absence of in-shock cooling. The ratio of the 
cooling-time to shock crossing time is of order 1.2. 



is that spatial location is not an accurate indication of the 
cooling time, as the latter is dependent upon the density. 
We would not necessarily expect the shock-tube simulations 
to have the same density profile or the same density jump 
at the shock front as the theoretical profile. However they 
do converge to the correct jump conditions far downstream. 

The sort of objects in numerical simulations where in- 
shock cooling will be important are those with a low ratio 
of cooling time to shock crossing time and which cool on a 
time-scale of order the dynamical time of the virialized halo. 
Objects which cool on a time-scale much longer than this are 
generally of little interest as they will not collapse sufficiently 
quick to form luminous objects. At the other extreme, if an 
object cools on a time-scale much less than the dynamical 
time the gas will be in a state of free-fall collapse; in-shock 
cooling will certainly be present in this case as the cooling 
is so rapid, but the correct conditions will be recovered on 
a time-scale less than the dynamical time. 

As an example, consider the collapse of an object with 
a cooling time equal to its collapse time, in a cosmological 
simulation with a dark-matter particle mass of 3 x 10 10 M o . 
Taking Q,b = 0.06, a galaxy with total mass of 1.7 x 1O 12 M0 
collapsing at redshift 2 would have a shock velocity at about 
350kms _1 . Taking the shock to be spread over 3 inter- 
particle spacings, we find the ratio of the cooling time to 
the shock crossing time to be 1.2. If the same object were to 
be simulated with a particle mass 1000 times lower, the ratio 
would increase to 12. The shock-tube simulations described 
in this section suggest that, for either large-scale cosmologi- 
cal simulations with a ratio of 1.2, or simulations of a single 
galaxy with a ratio of 8 the effect of in-shock cooling will 
reduce the maximum temperature reached by between 25 
to 50 percent and reduce the cooling time by a factor of 
approximately 2. 



5 SIMULATIONS WITH ELECTRON 
COOLING 

5.1 Theory 

The most important cooling mechanisms for a primordial 
gas with a temperature above 10 4 K result from recombina- 
tion of electrons with hydrogen ions and collisional excita- 
tion of neutral hydrogen by collisions with free electrons, the 
excited hydrogen atom then emits a photon which, should it 
not be re-absorbed, amounts to a reduction in the internal 
energy of the gas. The cooling terms, taken from Haiman, 
Thoul and Loeb (1996) are: 
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where T 5 is the temperature in units of 10 5 K. The rate of 
energy loss from the gas will be proportional to the num- 
ber densities of both neutral hydrogen and free electrons. 
The ionization level thus needs to be known and is gener- 
ally not equal to its equilibrium value for the situations we 
shall consider. In order to trace the ionization level we need 
to specify the ionization and recombination rates. For the 
purpose of this section we consider a gas composed entirely 
of hydrogen which is ionized through the reaction: 



H + e~ i — > H + + 2e~ 

This proceeds at a rate (Black 1978): 
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The gas recombines via the reaction: 
H + + e~ i — > H + hv 

which proceeds at a rate (Hutchins 1976): 
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5.2 Methodology 

For the electron-cooling runs, the gas was started with jump 
conditions for a perfectly isothermal shock. The motivation 
for this approach is that the cooling time is short, meaning 
that the perfectly isothermal conditions are reasonably close 
to the stable profile and should soon evolve to that state. 
The mass of the downstream particles was adjusted so that 
the correct density jump was obtained and the simulations 
were then allowed to run until the developing temperature, 
density and ionization profiles reach a stable state. Gas up- 
stream had an ionization level of 4 x 10~ 4 and a temperature 
of 12 000 K, which was the minimum temperature permitted 
in the simulations. Gas at or below 12 000 K did not have its 
ionization level updated and was not allowed to cool further 
as this would slow down the code too much. Instead, cooling 
below 12 000 K was estimated analytically. 

Run 1 is a fiducial run with units designed to represent 
an object at a redshift of 20 containing 10 6 particles with a 
total baryonic mass of 4 x 1O 8 M and a virial temperature 
of 78 000 K. Runs 2 and 3 have the same mass resolution as 
the fiducial run but represent objects with higher and lower 
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virial temperatures respectively. Runs 4 and 6 have the same 
virial temperature as the fiducial run but lower and higher 
mass resolution respectively. Runs 5 and 7 have a lower virial 
temperature than the fiducial but different mass resolutions. 
The mass resolution is represented by the number of parti- 
cles an object of mass 4 x 10 s M© would contain, although 
the simulations which have been performed have the same 
number of particles in every case. 

5.3 Results 

Table 1 below lists the results of seven shock-tube simula- 
tions using the cooling functions of equations [li] and |l9[ 
From left to right, the column headings of Table 1 are: run 
number; the number of particles an object of mass 4 x 1O 8 M0 
would contain, iV bj; the (analytic) virial temperature T a 
(the temperature the shocked gas should jump to); the ratio 
of the maximum temperature obtained in the simulations to 
the virial temperature, T 3 /T a ; the analytic ionization level 
for gas cooling from the virial temperature to 12 000 K, x a \ 
the ionization level of gas in the simulation that has cooled 
back to 12 000 K divided by the analytic level, x s /x a ; the ra- 
tio of the simulated to analytic cooling times from the peak 
of the shock, down to 12 000 K, r s /r a ; the ratio of the nu- 
merical to the analytic cooling times to reach 600 K, r a / /r a f, 
and the ratio of the numerical to the analytic H2 fractions 
at 600 K, H 2 , s /H 2 ,a. 

The analytic ionization level at 12 000 K, x a , assumes 
an adiabatic shock followed by rapid cooling and was calcu- 
lated by integrating the shock conservation equations 3, ^ 
and [tl] with the cooling functions of equations [l^ and fus- 
ing non-equilibrium chemistry. Cooling below 12 000 K was 
calculated allowing the gas to cool at a constant density to 
600K via the H2 cooling mechanism. A detailed treatment 
of the gas chemistry was used incorporating formation and 
destruction of H2, ionization and recombination, the cooling 
functions of equations [l^ and [[9] and the H2 cooling func- 
tion. The destruction and cooling rates of H2 was taken from 
Martin, Schwarz and Mandy (1996). The simulated cooling 
time from 12 000 K to 600 K, r s f, and H 2 fraction ti2 S use 
x 3 as the starting ionization level; T a f and H2 a use x a as the 
input ionization. 

5.4 Discussion 

5.4-1 Fiducial Run 

The effect of in-shock cooling is dramatic and striking in the 
fiducial run. The maximum temperature the shock reaches 
is only 0.27 of its theoretical value. This leads to an in- 
crease in the cooling time for the gas to cool back down to 
12 000 K by an astonishing factor of 236. There are a num- 
ber of contributing factors to this extreme figure. Firstly, 
although the ionization level at 12 000 K is 9 times lower in 
the shocktube than expected, this is a conservative estimate 
of the difference, as the factor is somewhat higher than this 
at other points on the cooling tail. This is because the lower 
jump temperature obtained results in a much slower rate of 
of ionization. Also, the density of the theoretical profile is 
somewhat higher at any given temperature (above 12 000 K) 
as the gas has reduced its temperature by a greater factor 
than gas at the same temperature in the shock tube. e.g. 



1 — 1 — 1 — ' — 1 — 1 — ' — ' — 1 — 1 — r 
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Figure 3. (a) Variation of Ionization Level with time for gas 
passing through the shock of the fiducial run (dashed line) and 
run 7 (solid line), (b) Variation of Ionization Level with time for 
gas passing through the shock of run 3 (dashed line) and run 5 
(solid line). 



for the fiducial run, at 22 000 K the theoretical profile has 
a density of 16 while gas at this temperature in the shock- 
tube has has a density of only 3.5. These effects combined 
produce the high factors of T s /r a seen in Table 1. The con- 
sequences of this are of extreme importance as simulations 
of such objects will be meaningless unless in-shock cooling 
can be removed, as the increase in the cooling times of such 
objects will slow the collapse. 

Also of importance is the ability of an object to cool to 
temperatures below 12 000 K when molecular hydrogen be- 
comes the dominant coolant. Molecular hydrogen is formed 
via the reactions 

H + e~ 1 — > H~ + hv 

rT +H1 — .Ha + e- 

Its production is thus proportional to the ionization level, 
the free electrons acting as a catalyst. The effect of in-shock 
cooling on the fiducial run was found to halve the H2 abun- 
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cl No. 




T a /K 


T 3 /T a 






T sf/ T af 


r s /T a 


H 2 , s /H 2 , a 


1 


10 6 


78 000 


0.27 


0.14 


0.11 


236 


2.14 


0.50 


2 


10 6 


168 000 


0.15 


0.58 


0.08 


512 


2.01 


0.51 


3 


10 6 


39 000 


0.45 


0.029 


0.24 


48 


1.74 


0.61 


4 


10 3 


78 000 


0.22 


0.14 


0.09 


941 


2.40 


0.45 


5 


10 3 


39 000 


0.44 


0.029 


0.16 


139 


2.10 


0.52 


6 


10 9 


78 000 


0.39 


0.14 


0.24 


42 


1.63 


0.64 


7 


10 9 


39 000 


0.81 


0.029 


0.41 


9 


1.40 


0.74 



Table 1. Steady state properties of shock tube simulations with cooling due to collisional excitation of free electrons with neutral 
hydrogen. 



dance and to increase the cooling time from 12 000 K to 
600 K by a factor of two. 

5.4-2 Effect of Virial Temperature 

Increasing the shock jump temperature, whilst retaining the 
same spatial and mass resolution results in a modest increase 
in the maximum temperature obtained at the shock front. 
This difference is sufficient to cause the ionization level of the 
cooled downstream gas, x B to be higher than that obtained 
in the fiducial run by a factor of almost three. However, the 
correct ionization level is still some way from being obtained 
in each of runs 1, 2 and 3. Comparing these runs shows that 
the fractional temperature reached increases with decreasing 
virial temperature, and the ionization level gets closer to its 
theoretical value. However, these trends do not appear to 
carry over to the cooling times and H2 abundances when 
the virial temperature increases from 78 000 K to 16 8000 K. 
This is due to the complex interaction between reaction rates 
and cooling times. 

5-4-3 Effect of resolution 

As one might expect, improving the mass resolution reduces 
the severity of any in-shock effects. For jump temperatures 
of 78 000 K and 39 000 K, as the mass resolution is improved 
there is a steady convergence towards the correct values of 
the H2 abundance, the cooling time to 600K, the cooling 
time from the virial temperature to 12 000 K and the jump 
temperature. Essentially the convergence of the first three 
is determined by the convergence towards the correct jump 
temperature, as it is this temperature that determines the 
extent to which the gas is re-ionized which in turn deter- 
mines the rate of H2 production and the cooling time. Fig- 
ure U shows the ionization level as a function of time for an 
analytical shock and for three shocks each with a virial tem- 
perature of 39 000 K, but with different mass resolutions. In 
each case we have defined time zero to be the time when 
gas reaches the maximum temperature. Negative time rep- 
resents gas which is being heated by the shock but has not 
yet reached the peak. These figures show the effect of lower- 
ing the resolution from the analytical profile, through runs 
7, 3 and 5. Note that the y-scale is the same in each panel 
but the x-scale is not. 

We can see from the graphs that as resolution is im- 
proved there is a convergence of the final ionization level 
and cooling time towards the theoretical profile. However, 
even if an object were to contain 10 9 particles the effects of 



in-shock cooling would still remain significant, most notice- 
ably for r s /r a . Such resolutions are at present not possible, 
and even if they were the quantities we have considered are 
yet to converge. We must conclude that the resolution of 
typical simulations being performed is hopelessly inadequate 
and if sensible results are to be obtained a way of removing 
in-shock cooling must be found. 



6 TEMPERATURE OSCILLATIONS OF 
UNSTABLE SHOCKS. 

6.1 Theory 

For a planar radiative shock with a cooling function A oc T a , 
Chevalier and Imamura (1982) showed that for a < 0.4 the 
shock is unstable against perturbations and will undergo os- 
cillations both spatially and in the jump temperature. The 
latter results from the fact that as the shock front moves, 
the relative velocity of the upstream gas in the frame of 
the shock varies. This in turn results in a different solution 
of the jump equations. The situation is similar for stable 
shocks also: if perturbed, they also undergo oscillations but 
with a decaying amplitude. For values of a between —1 and 
2 Chevalier and Imamura found the frequency of oscillation, 
lo, to vary between 0.26 and 0.31 with units of Ui„/x s where 
x 3 is the mean distance from the shock front to the point 
upstream where the gas has cooled to its upstream tem- 
perature, and Uin is the velocity of the upstream gas. The 
maximum displacement of the shock front from its mean 
position, x max , was also calculated and was found to be 
Xmax = t coo iv s , where t coo i is the steady state cooling time 
and v s the speed the shock front is moving at. 

6.2 Methodology 

The following simulations verify the instabilities found by 
Chevalier and Imamura (1982) for the case of a realistic 
cooling function composed of a number of power law terms, 
(some of which are stable). An analytic profile of the cooling 
tail was calculated as in Section 4.2. The gas upstream had 
a temperature of 50000K with a density jump of 100 for 
the downstream gas. This results in a temperature jump of 
a factor of 19.7 to 9.8 x 10 5 K. The cooling function used 
was that of Section 4, meaning that gas with a temperature 
above 10 5 K is in the unstable regime. The fact that the 
shock is unstable means that a steady state will never be 
reached so the simulation was ran until the nature of the 
oscillations became apparent. 
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Figure 4. Position of shock front as a function of time for an 
unstable radiative shock. 
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Figure 5. Oscillations of the jump temperature of an unstable 
shock as a function of time. 



6.3 Results 

The unstable nature of both the position and jump temper- 
ature of radiative shocks is well demonstrated upon consid- 
eration of Figures ^| and [| Although the cooling function 
used is not strictly a power law, but rather a combination of 
power laws, we find that both the magnitude of the spatial 
oscillations, and their period, are in reasonable agreement 
with the work of Chevalier and Imamura as presented above. 
The frequency of oscillations u was found to be u) — 0.31, 
where u) = (2ir/P) (x s /ui n ) and P is the observed period of 
oscillations. The majority of heated gas in the simulation has 
q ~ —0.8. For the case of a — —1, Chevalier and Imamura 
found the oscillation frequency to be to — 0.26. The agree- 
ment is reasonable, although interpolation of the Chevalier 
and Imamura values to intermediate values of ui is somewhat 
uncertain as there is no obvious linear trend in their results. 



The maximum displacement of the shock was found to be 
x m ax = 8.0 which agrees well with the theoretical value of 

%max — T.8. 

The effects of in-shock cooling on the simulations is un- 
certain but, given that we know it is present, it may account 
for the extreme oscillations down to an isothermal transition 
as seen in Figure 5. We note, however, that the period and 
spatial magnitude of the oscillations appears to be correct 
even when in-shock cooling is present. In the absence of in- 
shock cooling, any halo with a virial temperature of more 
than a few times 10 K (i.e. galaxies) will be subject to such 
an instability. 



7 CONCLUSION 

In this paper we have clearly demonstrated that in-shock 
cooling has important consequences for numerical simula- 
tions of galaxy formation, especially at high redshifts, where 
collisional excitation of hydrogen by free electrons, cooling 
by recombination and cooling by molecular hydrogen are 
the relevant coolants, but also for galaxies with higher virial 
temperatures collapsing at lower redshifts. 

The effect of in-shock cooling on a stable radiative shock 
was found to decrease the maximum temperature reached 
by 25 percent when the ratio of the cooling-time to shock 
crossing time was 8. This value increased to 50 percent for 
a ratio of 1.2. The result of this was that the profiles of the 
cooling tails were found to have different spatial extents from 
the theoretical profile. From this we conclude that although 
the effects of in-shock cooling are quite small, the collapse 
rate of such objects could be altered along with the spatial 
extent of any cooling tail. 

The effect on primordial objects of in-shock cooling was 
found to be extremely significant: the ionization level emerg- 
ing from the shock was significantly lower by a factor of 2.4 
in the best resolved case, and up to a factor of 12.1 in the 
worst case. The implications for the cooling time of the gas 
ranged from an increase of a factor of 9 up to nearly 1000, for 
the cooling times above 12 000 K. Below 12 000 K the result 
was less extreme, ranging from a factor of 1.4 to 2.4. This dif- 
ference is nevertheless important and will significantly slow 
the collapse of any object. The resulting molecular hydrogen 
abundance of the gas after it was allowed to cool to 600K 
was found in all cases to decrease from its analytic value by 
a factor of around 2. 

For the case of an unstable shock we find oscillations 
of the spatial location of the shock front and the maximum 
temperature reached. It is likely that the extreme oscillations 
of temperature observed (cooling to a virtually isothermal 
state) are fueled by in-shock cooling. For the resolutions 
currently achievable, an incorrect result will be obtained for 
simulations of galaxy formation. 

Finally, Figure g shows a successful attempt to remove 
in-shock cooling from the fiducial run. The method used to 
remove the effect was case-specific and was achieved by al- 
lowing the gas to cool only when its temperature is in excess 
of the adiabatic jump temperature or its density exceeds the 
adiabatic jump density. The upper curve in Figure ^| shows 
the gas density, the lower curve shows the temperature and 
the middle curve, with the double maximum, shows the di- 
mensionless ratio of the artificial pressure to the pressure, 
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Figure 6. Graph showing the properties of temperature, density 
and the ratio of the artificial pressure to the pressure for a shock 
where in-shock cooling has been removed. The temperature starts 
at 0.054 and is heated to 0.30 at position 56, this temperature 
being the correct jump temperature. The gas then cools extremely 
rapidly back down to its upstream value. The density starts at 
1 and increases as the gas shocks and the cools until it reaches 
its downstream value of 12. The ratio q/p has a double peaked 
structure, each peak being caused by shock heating and rapid 
cooling respectively. 
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q/p. The density and temperature are plotted in code units. 
The finite width of the shock can clearly be seen and the cor- 
rect jump conditions and ionization level both at the peak of 
the shock and downstream are recovered. The correct cool- 
ing tail down to a temperature of 12 000 K is also recovered. 

The importance of the ratio q/p is that it is a dimension- 
less indicator of when the gas is being shocked and might 
therefore be used as a method of removing in-shock cool- 
ing from other simulations, by means of preventing the gas 
from cooling when q/p exceeds a certain threshold value (i.e. 
when the gas is being shocked). However, this approach is 
unlikely to prove successful for two reasons. Firstly, consid- 
eration of the double maximum structure of q/p as seen in 
Figure ^| means than gas which is shocking does not nec- 
essarily have a value of q/p greater than gas which is not. 
Indeed, gas which is radiating rapidly (which is the case for 
gas whose density is increasing) has values of q/p as high as 
gas which is being shocked. Secondly, although q/p is dimen- 
sionless and in theory could produce a threshold value true 
for all scales of simulations, we would not expect the thresh- 
old value to be the same for shocks of different strengths. 
These two reasons rule out the use of q/p as a possible means 
of eliminating in-shock cooling. This was tested in a large- 
scale cosmological simulation. Different objects were found 
to have very different values of q/p and it was obvious that 
there was no single value that could be used to determine 
when gas was being shocked. 
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